{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "1260ea17-8840-46b7-a208-053c6ed797fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from sklearn.model_selection import train_test_split\n",
    "from nltk.corpus import stopwords\n",
    "from nltk.stem import WordNetLemmatizer\n",
    "import nltk\n",
    "import re\n",
    "import os\n",
    "import numpy as np\n",
    "from nltk import download\n",
    "from gensim.models import LdaMulticore, CoherenceModel\n",
    "from gensim.corpora.dictionary import Dictionary\n",
    "from nltk.tokenize import word_tokenize\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "5bcf9e89-6cd1-42d5-82ce-68db4e17e5dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('/zfs/disinfo/dcweekly_contrast/story_text/dc_weekly_topic_scores.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "86f15da0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove rows with month 2023-11\n",
    "df = df[df['month'] != '2023-11']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "9027a3f8-384d-440c-a702-cf119128fd86",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a binary indicator for whether we are in the AI-generated text period of time\n",
    "df['period'] = df.month.apply(lambda x: 1- int(x=='2023-06'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "c96a6ed5-5e6e-4fcd-b22d-78de199ef1d0",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[nltk_data] Downloading package stopwords to\n",
      "[nltk_data]     /scratch/cehrett/nltk_data...\n",
      "[nltk_data]   Package stopwords is already up-to-date!\n",
      "[nltk_data] Downloading package wordnet to\n",
      "[nltk_data]     /scratch/cehrett/nltk_data...\n",
      "[nltk_data]   Package wordnet is already up-to-date!\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Explicitly set the NLTK data directory\n",
    "nltk_data_dir = os.path.join('/', 'scratch', 'cehrett', 'nltk_data')\n",
    "os.environ['NLTK_DATA'] = nltk_data_dir\n",
    "\n",
    "# Ensure the directory is added to NLTK's search path\n",
    "nltk.data.path.append(nltk_data_dir)\n",
    "\n",
    "# Download the necessary packages\n",
    "nltk.download('stopwords', download_dir=nltk_data_dir)\n",
    "nltk.download('wordnet', download_dir=nltk_data_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "b93a9fab-776b-4793-a4fa-872d4e78b3b0",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[nltk_data] Downloading package stopwords to\n",
      "[nltk_data]     /scratch/cehrett/nltk_data...\n",
      "[nltk_data]   Package stopwords is already up-to-date!\n",
      "[nltk_data] Downloading package wordnet to\n",
      "[nltk_data]     /scratch/cehrett/nltk_data...\n",
      "[nltk_data]   Package wordnet is already up-to-date!\n",
      "[nltk_data] Downloading package punkt to /scratch/cehrett/nltk_data...\n",
      "[nltk_data]   Package punkt is already up-to-date!\n"
     ]
    }
   ],
   "source": [
    "# Download necessary NLTK data\n",
    "download('stopwords', download_dir=os.path.join('/', 'scratch', 'cehrett', 'nltk_data'))\n",
    "download('wordnet', download_dir=os.path.join('/', 'scratch', 'cehrett', 'nltk_data'))\n",
    "download('punkt', download_dir=os.path.join('/', 'scratch', 'cehrett', 'nltk_data'))\n",
    "\n",
    "# Function to clean and preprocess text\n",
    "def preprocess_text(text):\n",
    "    # Remove HTML tags\n",
    "    text = re.sub(r'<.*?>', '', text)\n",
    "    # Remove punctuation and numbers\n",
    "    text = re.sub(r'[^a-zA-Z\\s]', '', text)\n",
    "    # Lowercase\n",
    "    text = text.lower()\n",
    "    # Remove stopwords\n",
    "    stop_words = set(stopwords.words('english'))\n",
    "    text = ' '.join([word for word in text.split() if word not in stop_words])\n",
    "    # Lemmatization\n",
    "    lemmatizer = WordNetLemmatizer()\n",
    "    text = ' '.join([lemmatizer.lemmatize(word) for word in text.split()])\n",
    "    return text\n",
    "\n",
    "# Preprocess the text in the DataFrame\n",
    "df['processed_text'] = df['clean_text'].apply(preprocess_text)\n",
    "\n",
    "# Split the DataFrame into two categories based on the 'period' column\n",
    "df_period_0 = df[df['period'] == 0]\n",
    "df_period_1 = df[df['period'] == 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "de788165-987d-48dd-bf3c-fc4e45a0a176",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_3297322/2504636202.py:2: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_period_0['tokens'] = df_period_0['processed_text'].apply(word_tokenize)\n",
      "/tmp/ipykernel_3297322/2504636202.py:3: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  df_period_1['tokens'] = df_period_1['processed_text'].apply(word_tokenize)\n"
     ]
    }
   ],
   "source": [
    "# Tokenize the texts\n",
    "df_period_0['tokens'] = df_period_0['processed_text'].apply(word_tokenize)\n",
    "df_period_1['tokens'] = df_period_1['processed_text'].apply(word_tokenize)\n",
    "\n",
    "# Create Dictionaries and Corpora\n",
    "dictionary_0 = Dictionary(df_period_0['tokens'])\n",
    "corpus_0 = [dictionary_0.doc2bow(text) for text in df_period_0['tokens']]\n",
    "\n",
    "dictionary_1 = Dictionary(df_period_1['tokens'])\n",
    "corpus_1 = [dictionary_1.doc2bow(text) for text in df_period_1['tokens']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b0d3a5d3-c122-4476-88c9-e893278729d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_coherence_values(dictionary, corpus, texts, start, limit, step):\n",
    "    \"\"\"\n",
    "    Computes coherence scores for LDA models with varying number of topics.\n",
    "    \n",
    "    This function iterates over a specified range of topic numbers, fitting an LDA model for each count,\n",
    "    and calculates the coherence score to evaluate the model's quality. It uses the LdaMulticore class\n",
    "    for efficient parallel computation. Key hyperparameters are chosen to balance model complexity,\n",
    "    convergence speed, and coherence quality:\n",
    "    \n",
    "    - `num_topics`: Varies in steps from `start` to `limit` to explore the optimal number of topics.\n",
    "    - `random_state`: Ensures reproducibility of results.\n",
    "    - `chunksize`: Number of documents to process at a time in the training algorithm. A balance between\n",
    "      memory consumption and computational efficiency.\n",
    "    - `passes`: The number of passes through the corpus during training, for sufficient model convergence.\n",
    "    - `alpha`: Set to 'asymmetric' to account for the intuition that some topics are more probable than others.\n",
    "    - `eta`: Left as `None` to use the default symmetric prior, assuming no prior knowledge about topic-word distribution.\n",
    "    - `decay` and `offset`: Control the learning rate, helping to stabilize early iterations of the model training.\n",
    "    - `eval_every`: Determines how often to calculate log perplexity; set to optimize the trade-off between\n",
    "      training speed and model evaluation.\n",
    "    - `iterations`: Max number of iterations over each document; higher for better model precision.\n",
    "    - `gamma_threshold`: Convergence threshold for gamma, the document-topic density.\n",
    "    - `minimum_probability`: Filters out topics with a probability below this threshold in the topics-per-document distribution.\n",
    "    - `minimum_phi_value`: Lower bound on term probabilities, filtering out terms.\n",
    "    - `per_word_topics`: If True, computes a list of most likely topics for each word, alongside its phi value times feature length.\n",
    "    - `workers`: Utilizes 7 of 8 available cores, leaving one free to prevent resource saturation.\n",
    "    \n",
    "    Parameters:\n",
    "        dictionary (gensim.corpora.dictionary.Dictionary): The dictionary mapping of id->word.\n",
    "        corpus (list of list of (int, float)): Corpus represented as a list of document vectors or sparse matrix.\n",
    "        texts (list of list of str): Text data for coherence calculation; should match `corpus`.\n",
    "        start (int): The starting number of topics.\n",
    "        limit (int): The max number of topics to model.\n",
    "        step (int): The step size to iterate through the number of topics.\n",
    "    \n",
    "    Returns:\n",
    "        model_list (list of gensim.models.LdaMulticore): List of LDA multicore models for each number of topics.\n",
    "        coherence_values (list of float): Coherence values corresponding to each model. \n",
    "        The `c_v` coherence measure is chosen for its balance between interpretability and informativeness. \n",
    "        Unlike other coherence measures that may rely solely on document co-occurrence data (`u_mass`) \n",
    "        or on similarity between top words (`c_uci`, `c_npmi`), `c_v` combines a variety of data sources, \n",
    "        including word co-occurrence and document-level information, to provide a robust assessment of topic quality. \n",
    "        This measure has been shown to correlate well with human judgments of topic coherence, \n",
    "        making it a suitable choice for evaluating the meaningfulness and distinctiveness of the topics discovered by LDA models. \n",
    "        Additionally, `c_v` allows for the incorporation of sliding window techniques, \n",
    "        which can capture more nuanced semantic relationships between words, further enhancing the evaluation of topic coherence.\n",
    "\n",
    "    \"\"\"\n",
    "    coherence_values = []\n",
    "    model_list = []\n",
    "    for num_topics in range(start, limit, step):\n",
    "        print(f'Beginning LDA with {num_topics} topics')\n",
    "        model = LdaMulticore(corpus=corpus,\n",
    "                             id2word=dictionary,\n",
    "                             num_topics=num_topics,\n",
    "                             random_state=100,\n",
    "                             chunksize=100,\n",
    "                             passes=10,\n",
    "                             alpha='asymmetric',\n",
    "                             eta=None,\n",
    "                             decay=0.5,\n",
    "                             offset=64,\n",
    "                             eval_every=10,\n",
    "                             iterations=400,\n",
    "                             gamma_threshold=0.001,\n",
    "                             minimum_probability=0.01,\n",
    "                             minimum_phi_value=0.01,\n",
    "                             per_word_topics=True)\n",
    "        model_list.append(model)\n",
    "        print('LDA model fitted. Beginning coherence calculation')\n",
    "        coherencemodel = CoherenceModel(model=model, texts=texts, dictionary=dictionary, coherence='c_v')\n",
    "        coherence_values.append(coherencemodel.get_coherence())\n",
    "\n",
    "    return model_list, coherence_values\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "808621a9-2995-4f44-b792-09ce08c32b1b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Beginning LDA with 4 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 8 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 12 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 16 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 20 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 24 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 28 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 32 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 36 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 40 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 44 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 48 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 52 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 56 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 60 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 64 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 68 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 72 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 76 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 80 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 84 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 88 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 92 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 96 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 100 topics\n",
      "LDA model fitted. Beginning coherence calculation\n"
     ]
    }
   ],
   "source": [
    "# Compute coherence values for each set of topics\n",
    "model_list_0, coherence_values_0 = compute_coherence_values(dictionary=dictionary_0, corpus=corpus_0, texts=df_period_0['tokens'], start=4, limit=101, step=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "4eef5672-36a7-45a5-9461-070501b55041",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Beginning LDA with 4 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 8 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 12 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 16 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 20 topics\n",
      "LDA model fitted. Beginning coherence calculation\n",
      "Beginning LDA with 24 topics\n",
      "LDA model fitted. Beginning coherence calculation\n"
     ]
    }
   ],
   "source": [
    "# Compute coherence values for each set of topics\n",
    "model_list_1, coherence_values_1 = compute_coherence_values(dictionary=dictionary_1, corpus=corpus_1, texts=df_period_1['tokens'], start=4, limit=101, step=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0903ed50-406d-4417-b7b6-0cb8660eb397",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = {\n",
    "    'set_0': {\n",
    "        'models': model_list_0,\n",
    "        'coherence_values': coherence_values_0\n",
    "    },\n",
    "    'set_1': {\n",
    "        'models': model_list_1,\n",
    "        'coherence_values': coherence_values_1\n",
    "    }\n",
    "}\n",
    "\n",
    "import pickle\n",
    "\n",
    "# Save the packaged object\n",
    "with open('/zfs/disinfo/dcweekly_contrast/story_text/lda_results.pkl', 'wb') as f:\n",
    "    pickle.dump(results, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "81a51eeb-aa3e-4d54-b364-6c3b3c58b006",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "# Load the packaged object\n",
    "with open('/zfs/disinfo/dcweekly_contrast/story_text/lda_results.pkl', 'rb') as f:\n",
    "    loaded_results = pickle.load(f)\n",
    "\n",
    "# Accessing loaded results\n",
    "model_list_0 = loaded_results['set_0']['models']\n",
    "coherence_values_0 = loaded_results['set_0']['coherence_values']\n",
    "model_list_1 = loaded_results['set_1']['models']\n",
    "coherence_values_1 = loaded_results['set_1']['coherence_values']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "c6c554b7-d50c-4adc-aeed-f213672786a2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwwAAAEiCAYAAABURlUUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbhklEQVR4nO3dd1xTV/8H8E/YG3GxVbQu3HWLiLio21r31ra2jtZda5fWR0Vt66N1W+sozgd3XXUAihsHbq17gVIn1IESzu+P80skJoEEAmF83q9XXpB7T26+uYSTfO9ZCiGEABERERERkQ4W5g6AiIiIiIhyLyYMRERERESkFxMGIiIiIiLSiwkDERERERHpxYSBiIiIiIj0YsJARERERER6MWEgIiIiIiK9mDAQEREREZFeTBiIiIiIiEgvJgxE+dyZM2fQv39/+Pn5wc7ODk5OTnj//fcxffp0PH782OjjNW7cGJUrV86GSPOPR48eYdy4cfD394ejoyNcXV1RoUIF9O7dG2fOnDF3eNli7dq1qFSpEuzt7aFQKBAbG5ttzxUVFQWFQoF169alW06hUKhvlpaWcHNzQ7Vq1fDZZ5/hyJEj6T72119/hUKhMPq93rhxYygUCpQuXRpCCK39+/fvV8e0bNkyo46dnmXLlkGhUODmzZtGP3bChAlQKBQmi4WI8h8mDET52G+//YaaNWsiJiYGY8aMwc6dO7Fx40Z07twZCxYswMcff2zuEPOdf//9F/Xq1cOyZcvwySefYMuWLVi5ciUGDhyIGzduZOsXaXP5559/0Lt3b5QpUwY7d+7E4cOHUa5cOXOHBQDo1KkTDh8+jAMHDmDNmjXo06cPjhw5gvr162PYsGF6H7dkyRIAwPnz53H06FGjntPZ2Rk3btxARESEzuO6uLgY9yKIiMzMytwBEFH2OHz4MAYNGoTmzZtj06ZNsLW1Ve9r3rw5Ro0ahZ07d5oxQv2EEHj16hXs7e3NHYrRwsPDcfXqVURERCA4OFhj38iRI5Gamppjsbx58wYKhQJWVtlb1f/999948+YNevXqhaCgIJMc88WLF3BwcMjycdzd3VGvXj31/ZCQEAwfPhwDBw7Er7/+igoVKmDQoEEajzl+/DhOnz6N1q1bY9u2bfj9999Rt25dg5+zRIkScHZ2xpIlS9C0aVP19qSkJISHh6Nnz5747bffsvzaiIhyClsYiPKpKVOmQKFQYNGiRRrJgoqNjQ3atWunvp+amorp06ejQoUKsLW1RfHixdGnTx/cvXtX5/FjYmIQGBgIBwcHlC5dGlOnTtX6MpyYmIjRo0fDz88PNjY28Pb2xvDhw/H8+XONcgqFAkOHDsWCBQtQsWJF2NraYvny5QCAK1euoEePHihevDhsbW1RsWJFzJ07V+Pxqi4qq1evxrfffgsvLy+4uLigWbNmuHz5slbsO3fuRNOmTeHq6goHBwdUrFgRoaGhGmWOHz+Odu3aoXDhwrCzs0ONGjXwv//9L50zLj169AgA4OnpqXO/hYVmtXvp0iV0794d7u7usLW1RYkSJdCnTx8kJyery5w7dw7t27eHm5sb7OzsUL16dfX5efcchIWFYdSoUfD29oatrS2uXr0KANizZw+aNm0KFxcXODg4ICAgAHv37tU4xj///IOBAwfC19cXtra2KFasGAICArBnzx69r7dfv35o2LAhAKBr165QKBRo3Lixev+WLVtQv359ODg4wNnZGc2bN8fhw4c1jqHqEnPy5El06tQJbm5uKFOmjN7nzCpLS0vMmTMHRYsWxU8//aS1//fffwcATJ06FQ0aNMCaNWvw4sULo55jwIAB2LBhA54+faretmbNGgBAt27ddD7mwIEDaNq0KZydneHg4IAGDRpg27ZtWuWOHDmCgIAA2NnZwcvLC+PGjcObN290HnPt2rWoX78+HB0d4eTkhJCQEJw6dSrD+CMiItC4cWMUKVIE9vb2KFGiBD766COjzwMR5ROCiPKdlJQU4eDgIOrWrWvwYwYOHCgAiKFDh4qdO3eKBQsWiGLFiglfX1/xzz//qMsFBQWJIkWKiLJly4oFCxaI3bt3i8GDBwsAYvny5epyz58/F9WrVxdFixYVM2bMEHv27BGzZs0Srq6uokmTJiI1NVVdFoDw9vYWVatWFatWrRIRERHi3Llz4vz588LV1VVUqVJF/PHHH2LXrl1i1KhRwsLCQkyYMEH9+MjISAFAlCpVSvTs2VNs27ZNrF69WpQoUUKULVtWpKSkqMsuXrxYKBQK0bhxY7Fq1SqxZ88eMW/ePDF48GB1mYiICGFjYyMCAwPF2rVrxc6dO0W/fv0EALF06dJ0z+OBAwcEAFG7dm2xceNG8fDhQ71lY2NjhZOTkyhVqpRYsGCB2Lt3r1ixYoXo0qWLSExMFEIIcenSJeHs7CzKlCkj/vjjD7Ft2zbRvXt3AUBMmzZN6xx4e3uLTp06iS1btoitW7eKR48eibCwMKFQKESHDh3Ehg0bxJ9//inatGkjLC0txZ49e9THCAkJEcWKFROLFi0SUVFRYtOmTeKHH34Qa9as0fsarl69KubOnSsAiClTpojDhw+L8+fPCyGEWLlypQAgWrRoITZt2iTWrl0ratasKWxsbER0dLT6GOPHjxcARMmSJcXYsWPF7t27xaZNm/Q+p+q1hoeHp/u3ACCGDBmid3+3bt0EAHHnzh31thcvXghXV1dRu3ZtIYR8vwAQy5YtS/e5VIKCgkSlSpVEYmKicHR0FPPmzVPvq1u3rujTp4+IiYnRei9FRUUJa2trUbNmTbF27VqxadMm0aJFC6FQKDTO//nz54WDg4Pw9/cXq1evFps3bxYhISGiRIkSAoC4ceOGuuzkyZOFQqEQAwYMEFu3bhUbNmwQ9evXF46Ojuq/kRBvz7/KjRs3hJ2dnWjevLnYtGmTiIqKEitXrhS9e/cWT548Meg8EFH+woSBKB+6f/++ACC6detmUPmLFy8KABpfmoUQ4ujRowKA+Oabb9TbgoKCBABx9OhRjbL+/v4iJCREfT80NFRYWFiImJgYjXLr1q0TAMT27dvV2wAIV1dX8fjxY42yISEhwsfHRzx79kxj+9ChQ4WdnZ26vOoLZKtWrTTK/e9//xMAxOHDh4UQQiQlJQkXFxfRsGFDjYTlXRUqVBA1atQQb9680djepk0b4enpKZRKpd7HCiHExIkThY2NjQAgAAg/Pz/x+eefi9OnT2uUa9KkiShUqJBISEjQe6xu3boJW1tbcfv2bY3tLVu2FA4ODuLp06ca56BRo0Ya5Z4/fy4KFy4s2rZtq7FdqVSKatWqiTp16qi3OTk5ieHDh6f72nTR9QVeqVQKLy8vUaVKFY3zlZSUJIoXLy4aNGig3qb6wvrDDz9k+vl0yShhGDt2rNZ7+Y8//hAAxIIFC9TxOjk5icDAQINiUyUMQgjRt29fUatWLSGE/KIPQERFRelMGOrVqyeKFy8ukpKS1NtSUlJE5cqVhY+Pj/r92rVrV2Fvby/u37+vUa5ChQoaCcPt27eFlZWV+OKLLzTiS0pKEh4eHqJLly7qbe8mDKr/0djYWINeMxHlf+ySRESIjIwEILuXpFWnTh1UrFhRq+uKh4cH6tSpo7GtatWquHXrlvr+1q1bUblyZVSvXh0pKSnqW0hICBQKBaKiojQe36RJE7i5uanvv3r1Cnv37sWHH34IBwcHjWO0atUKr1690prpJm0XK1VMANRxHTp0CImJiRg8eLDeWWGuXr2KS5cuoWfPngCg9bzx8fE6uzml9f333+P27dtYsmQJPvvsMzg5OWHBggWoWbMmVq9eDUD20d+3bx+6dOmCYsWK6T1WREQEmjZtCl9fX43t/fr1w4sXL7S693z00Uca9w8dOoTHjx+jb9++Gq8lNTUVH3zwAWJiYtRdxOrUqYNly5Zh0qRJOHLkiN5uLoa4fPky4uLi0Lt3b41uWE5OTvjoo49w5MgRre4t78ae3YSOWYx+//132Nvbq7sNOTk5oXPnzoiOjsaVK1eMOv6AAQNw/PhxnD17Fr///jvKlCmDRo0aaZV7/vw5jh49ik6dOsHJyUm93dLSEr1798bdu3fV77nIyEg0bdoU7u7uGuW6du2qccy//voLKSkp6NOnj8bf3c7ODkFBQVr/f2lVr14dNjY2GDhwIJYvX47r168b9bqJKP9hwkCUDxUtWhQODg64ceOGQeXT63fv5eWl3q9SpEgRrXK2trZ4+fKl+v6DBw9w5swZWFtba9ycnZ0hhMDDhw81Hv/ucz969AgpKSmYPXu21jFatWoFAFrHeDcu1dgNVVz//PMPAMDHx0fPmZBxA8Do0aO1nnfw4ME6n1cXd3d39O/fHwsWLMCZM2ewb98+2NjYqGfmefLkCZRKZbqxqM6Dvr+Lan9a75ZVvZ5OnTppvZ5p06ZBCKGeXnft2rXo27cvFi9ejPr166Nw4cLo06cP7t+/n+Hr1RW3rnhUsaempuLJkyfpxp7dVImk6lxevXoV+/fvR+vWrSGEwNOnT/H06VN06tQJwNuZkwzVqFEjlC1bFgsXLkRYWBgGDBigM1F98uQJhBAG/Z0fPXoEDw8PrXLvblP93WvXrq31d1+7dm267+EyZcpgz549KF68OIYMGYIyZcqgTJkymDVrluEvnojyFc6SRJQPWVpaomnTptixYwfu3r2b4ZdS1Rft+Ph4rbJxcXEoWrSo0TEULVoU9vb2er9kvXvMd79Iubm5qa+wDhkyROcx/Pz8jIpJdSVf30DutHGNGzcOHTt21FmmfPnyRj0vIL88tmjRAps2bUJCQgIKFy4MS0vLdGMB5N8mPj5ea3tcXJxGvCrvnkfV/tmzZ2vMFpSW6mp10aJFMXPmTMycORO3b9/Gli1b8PXXXyMhIcHoGbXSvqd0xW5hYaHRoqQr9uz08uVL7NmzB2XKlFG/55csWQIhBNatW6dzjYfly5dj0qRJsLS0NPh5+vfvj++++w4KhQJ9+/bVWcbNzQ0WFhYG/Z2LFCmiM4F7d5uq/Lp161CyZEmD41UJDAxEYGAglEoljh8/jtmzZ2P48OFwd3fXO2ibiPIvJgxE+dS4ceOwfft2fPrpp9i8eTNsbGw09r958wY7d+5E27Zt0aRJEwDAihUrULt2bXWZmJgYXLx4Ed9++63Rz9+mTRtMmTIFRYoUMfqLPQA4ODggODgYp06dQtWqVbXiz4wGDRrA1dUVCxYsQLdu3XR+QS1fvjzKli2L06dPY8qUKUY/x4MHD1CsWDGt2ZCUSiWuXLkCBwcHFCpUCDY2NggKCkJ4eDgmT56sNylr2rQpNm7ciLi4OPXVZgD4448/4ODgoDcJUAkICEChQoVw4cIFDB061ODXUaJECQwdOhR79+7FwYMHDX6cSvny5eHt7Y1Vq1Zh9OjR6nP9/PlzrF+/Xj1zkjkolUoMHToUjx49Us+OpVQqsXz5cpQpUwaLFy/WeszWrVvxyy+/YMeOHWjTpo3Bz9W3b18cPXoUFStWhLe3t84yjo6OqFu3LjZs2ICff/5ZPZ1wamoqVqxYAR8fH/W6FsHBwdiyZQsePHigTvSUSiXWrl2rccyQkBBYWVnh2rVrWerqZWlpibp166JChQpYuXIlTp48yYSBqABiwkCUT9WvXx/z58/H4MGDUbNmTQwaNAiVKlXCmzdvcOrUKSxatAiVK1dG27ZtUb58eQwcOBCzZ8+GhYUFWrZsiZs3b+L777+Hr68vRowYYfTzDx8+HOvXr0ejRo0wYsQIVK1aFampqbh9+zZ27dqFUaNGZTi3/axZs9CwYUMEBgZi0KBBKFWqFJKSknD16lX8+eefOhfGSo+TkxN++eUXfPLJJ2jWrBk+/fRTuLu74+rVqzh9+jTmzJkDAFi4cCFatmyJkJAQ9OvXD97e3nj8+DEuXryIkydPIjw8XO9zhIWFYeHChejRowdq164NV1dX3L17F4sXL8b58+fxww8/qJOfGTNmoGHDhqhbty6+/vprvPfee3jw4AG2bNmChQsXwtnZGePHj8fWrVsRHByMH374AYULF8bKlSuxbds2TJ8+Ha6urhm+5tmzZ6Nv3754/PgxOnXqhOLFi+Off/7B6dOn8c8//2D+/Pl49uwZgoOD0aNHD1SoUAHOzs6IiYnBzp079ba0pMfCwgLTp09Hz5490aZNG3z22WdITk7GTz/9hKdPn2Lq1KlGH/Nd+lZrDgoKUrcmPXjwAEeOHIEQAklJSTh37hz++OMPnD59GiNGjMCnn34KANixYwfi4uIwbdo0jWlhVSpXrow5c+bg999/Nyph8PLywqZNmzIsFxoaiubNmyM4OBijR4+GjY0N5s2bh3PnzmH16tXqhOu7777Dli1b0KRJE/zwww9wcHDA3LlztaYqLlWqFCZOnIhvv/0W169fxwcffAA3Nzc8ePAAx44dg6OjI3788UedsSxYsAARERFo3bo1SpQogVevXqlbCps1a2bwayeifMSMA66JKAfExsaKvn37ihIlSggbGxvh6OgoatSoIX744QeN2XmUSqWYNm2aKFeunLC2thZFixYVvXr10phyUgjNWWDS6tu3ryhZsqTGtn///Vd89913onz58sLGxkY9ReqIESM0ZnlBOrPZ3LhxQwwYMEB4e3sLa2trUaxYMdGgQQMxadIkdRl9s+bcuHFD51So27dvF0FBQcLR0VE9RWXaKUqFEOL06dOiS5cuonjx4sLa2lp4eHiIJk2aqGfP0efChQti1KhRolatWqJYsWLCyspKuLm5iaCgIBEWFqazfOfOnUWRIkWEjY2NKFGihOjXr5949eqVuszZs2dF27Zthaurq7CxsRHVqlXTek0ZzRy0b98+0bp1a1G4cGFhbW0tvL29RevWrdXlX716JT7//HNRtWpV4eLiIuzt7UX58uXF+PHjxfPnz9N9zek996ZNm0TdunWFnZ2dcHR0FE2bNhUHDx7UKKOapSft9L2GPJ++W2RkpBBCaGyzsLAQLi4uokqVKmLgwIHqmbNUOnToIGxsbDKcscrKykrjvfsuff8faemaJUkIIaKjo0WTJk2Eo6OjsLe3F/Xq1RN//vmn1uMPHjwo6tWrJ2xtbYWHh4cYM2aMWLRokda0qkLI8x8cHCxcXFyEra2tKFmypOjUqZPGdLrvzpJ0+PBh8eGHH4qSJUsKW1tbUaRIEREUFCS2bNmS7usiovxLIYSOaSKIiIiIiIjAWZKIiIiIiCgdTBiIiIiIiEgvJgxERERERKQXEwYiIiIiItKLCQMREREREenFhIGIiIiIiPRiwkBERERERHoxYSAiIiIiIr2YMBARERERkV5MGIiIiIiISC8mDEREREREpBcTBiIiIiIi0osJAxERERER6cWEgYiIiIiI9GLCQEREREREejFhICIiIiIivZgwEBERERGRXkwYiIiIiIhILyYMRERERESkFxMGIiIiIiLSiwkDERERERHpxYSBiIiIiIj0YsJARERERER6MWEgIiIiIiK9mDAQEREREZFeTBiIiIiIiEgvJgxERERERKQXEwYiIiIiItKLCQMREREREenFhIGIiIiIiPRiwkBERERERHoxYSAiIiIiIr2YMBARERERkV5MGIiIiIiISC8mDEREREREpBcTBiIiIiIi0osJAxERERER6cWEgYiIiIiI9GLCQEREREREejFhICIiIiIivZgwEBERERGRXkwYiExs2bJlUCgUOH78uLlDISIqUI4cOYLOnTvD09MTNjY28PDwQKdOnXD48GGjj3XhwgVMmDABN2/e1FvmzJkzUCgUOHXqFABAoVBAoVBg6tSpWmX52UB5GRMGIiIiyvNmz56NgIAA3L17F9OnT8eePXvw888/4969e2jYsCHmzJlj1PEuXLiAH3/8Md2EYf369fDz80ONGjU0tk+dOhWPHz/OzMsgypWYMBAREVGedvDgQQwfPhytWrVCdHQ0evfujUaNGqFXr16Ijo5Gq1atMGzYMBw8eNCkz7tu3Tp89NFHGtuaNWuG58+fY/LkySZ9LiJzYsJAlM0aN26Mxo0ba23v168fSpUqpb5/8+ZNKBQK/Pzzz5gxYwb8/Pzg5OSE+vXr48iRI1qPP378ONq1a4fChQvDzs4ONWrUwP/+979sfCVERLlTaGgoFAoF5s+fDysrK419VlZWmDdvnlZXoUuXLqF79+5wd3eHra0tSpQogT59+iA5ORnLli1D586dAQDBwcHqrkbLli3TePyFCxe0Eoby5cvj448/xty5c3Hr1q0MY9+yZQvq168PBwcHODs7o3nz5hpdqDZt2gSFQoG9e/dqPXb+/PlQKBQ4c+aMQeeJKLOYMBDlMnPnzsXu3bsxc+ZMrFy5Es+fP0erVq3w7NkzdZnIyEgEBATg6dOnWLBgATZv3ozq1auja9euGh9oRET5nVKpRGRkJGrVqgUfHx+dZXx9fVGzZk1ERERAqVTi9OnTqF27No4cOYKJEydix44dCA0NRXJyMl6/fo3WrVtjypQpAGSdfPjwYRw+fBitW7dWH3P9+vXw9vZG3bp1tZ5vwoQJsLS0xPfff59u7KtWrUL79u3h4uKC1atX4/fff8eTJ0/QuHFjHDhwAADQpk0bFC9eHEuXLtV6/LJly/D++++jatWqBp8vosywyrgIEeUkZ2dnbN26FZaWlgAALy8v1KlTBzt27EC3bt0AAIMHD0alSpUQERGhvpoWEhKChw8f4ptvvkGfPn1gYcHrAUSU/z18+BAvXryAn59fuuX8/Pxw7NgxPHr0CCNHjoSVlRWOHTuGYsWKqcv07NkTgKyHy5YtCwDw9/dHvXr1tI63bt06dOzYEQqFQmufh4cHRowYgdDQUIwePVrnF/rU1FSMGTMGVapUwY4dO9R1dqtWrVCmTBmMHTsWBw8ehJWVFXr16oX58+fj2bNncHV1BQBcvHgRx44dw+zZsw08U0SZx28URLlM69at1ckCAPUHjapp++rVq7h06ZL6gy0lJUV9a9WqFeLj43H58uWcD5yIKBcTQgAAXr58iX379qFLly4ayYIxrl+/jtjYWK3uSGl99dVXKFy4MMaOHatz/+XLlxEXF4fevXtrXOBxcnLCRx99hCNHjuDFixcAgAEDBuDly5dYu3atutzSpUtha2uLHj16ZOo1EBmDCQNRLlOkSBGN+7a2tgDkhxwAPHjwAAAwevRoWFtba9wGDx4MQF5xIyIqCIoWLQoHBwfcuHEj3XI3b96Eg4MDrKysoFQq9XZfMsS6detQvHhxNGzYUG8ZFxcXfPfdd9i5cyciIyO19j969AgA4OnpqbXPy8sLqampePLkCQCgUqVKqF27trpbklKpxIoVK9C+fXsULlw406+DyFBMGIiymZ2dHZKTk7W2Z/ZLfdGiRQEA48aNQ0xMjM5b9erVsxIyEVGeYWlpieDgYBw/fhx3797VWebu3bs4ceIEmjRpgsKFC8PS0lJvWUOsX78eHTp00GgN1mXQoEHw8/PD2LFj1S0cKqqLQ/Hx8VqPi4uLg4WFBdzc3NTb+vfvjyNHjuDixYvYuXMn4uPj0b9//0y/BiJjMGEgymalSpXC33//rZE0PHr0CIcOHcrU8cqXL4+yZcvi9OnTqFWrls6bs7OzqcInIsr1xo0bByEEBg8eDKVSqbFPqVRi0KBBEEJg3LhxsLe3R1BQEMLDw9O9cPNu667KnTt3EBMTk253JBUbGxtMmjQJMTExCA8P19hXvnx5eHt7Y9WqVRrJxPPnz7F+/Xr1zEkq3bt3h52dHZYtW4Zly5bB29sbLVq0yDAGIlNgwkCUTVQD4Xr37o3Hjx+jV69e2LVrF1avXo1mzZrBxcUl08deuHAh9u7di5CQEKxevRr79+/Hpk2bEBoaqp4KkIiooAgICMDMmTOxbds2NGzYECtXrkR0dDRWrlyJwMBAbN++HTNnzkSDBg0AADNmzMCbN29Qt25d/Pbbb4iMjMSaNWvQo0cPJCUlAQAqV64MAFi0aBEOHDiA48eP49GjR1i/fj0KFSqE4OBgg2Lr3r07atSogR07dmhst7CwwPTp0xEbG4s2bdpgy5YtCA8PR3BwMJ4+faq1WnShQoXw4YcfYtmyZdiyZQv69u3LyS0ox/CdRmRiqkFqqqtTAQEBWL58Oc6fP4/27dtj0qRJGDdunM61GQwVHByMY8eOoVChQhg+fDiaNWuGQYMGYc+ePWjWrJkpXgYRUZ7yxRdf4ODBg/Dx8cGoUaPQpEkTjBw5Ep6enjhw4AC++OILddlq1arh2LFjqFmzJsaNG4cPPvgAY8eOha2tLWxsbADIWZVmzpyJ06dPo3Hjxqhduzb+/PNPrF+/Hu3atYO1tbVBcSkUCkybNk3nvh49emDTpk149OgRunbtiv79+8PFxQWRkZE6x0f0798fCQkJeP36Nfr162f8SSLKJIV4t1MdEWXJsGHDMGfOHDx9+pRdg4iI8pH79+/D29sbmzZtQtu2bc0dDlGO4ToMRCZy4sQJxMTEYMmSJWjXrh2TBSKifMbDw0NrjARRQcAWBiIT8fPzw7Nnz9CyZUv8+uuvWtOjEhEREeVFTBiIiIiIiEgvDnomIiIiIiK9mDAQEREREZFeTBiIiIiIiEivAjdLUmpqKuLi4uDs7KxeWIuIiAwnhEBSUhK8vLyyfeEo1tlERJlnqvq6wCUMcXFx8PX1NXcYRER53p07d+Dj45Otz8E6m4go67JaXxe4hEE1N/6dO3fg4uJi5miIiPKexMRE+Pr65shaI6yziYgyz1T1dYFLGFRN2i4uLvzwISLKgpzoIpSlOlupBKKjgfh4wNMTCAwELC2zIUoiotwtq/V1gUsYiIioANiwARg2DLh79+02Hx9g1iygY0fzxUVElAdxliQiIspfNmwAOnXSTBYA4N49uX3DBvPERUSURzFhICKi/EOplC0LQmjvU20bPlyWIyIigzBhICKi/CM6WrtlIS0hgDt3ZDkiIjIIEwYiIso/4uNNW46IiJgwEBFRPuLpadpyRETEhIGIiPKRwEA5G5K+KQQVCsDXV5YjIiKDMGEgIqL8w9JSTp0K6E8aZs7kegxEREZgwkBERPlLx47AunWAt7f2vvnzuQ4DEZGRmDAQEVH+07EjcPMmEBkJrFwJVK8ut+/fb86oiIjyJCYMRESUP1laAo0bAz16AEuWyC5Kq1YBx46ZOzIiojyFCQMREeV/NWoAffrI30eO1L2wGxER6cSEgYiICobJkwF7e+DgQWDDBnNHQ0SUZzBhICKigsHbGxgzRv4+diyQnGzeeIiI8ggmDEREVHCMGQN4eADXrgFz55o7GiKiPIEJAxERFRxOTsCkSfL3//wHePTIvPEQEeUBmUoYwsLCEBAQAC8vL9y6dQsAMHPmTGzevNmkwREREZlcv35A1arA06cyaSAionQZnTDMnz8fI0eORKtWrfD06VMolUoAQKFChTBz5kxTx0dERGRalpbAL7/I3+fOBf7+27zxEBHlckYnDLNnz8Zvv/2Gb7/9FpaWlurttWrVwtmzZ00aHBERUbZo1gxo1QpISZEDoImISC+jE4YbN26gRo0aWtttbW3x/PlzkwRFRESU7X76SbY2bNoE7Ntn7miIiHItoxMGPz8/xMbGam3fsWMH/P39TRETERFR9vP3BwYOlL+PGgWkppo3HiKiXMrK2AeMGTMGQ4YMwatXryCEwLFjx7B69WqEhoZi8eLF2REjERFR9pgwAVixAjhxAli1CujVy9wRERHlOkYnDP3790dKSgq++uorvHjxAj169IC3tzdmzZqFbt26ZUeMRERE2aN4ceCbb4Bx4+StY0fAwcHcURER5SpGdUlKSUnB8uXL0bZtW9y6dQsJCQm4f/8+7ty5g48//ji7YiQiIso+w4cDJUoAd+8C//2vuaMhIsp1jEoYrKysMGjQICQnJwMAihYtiuLFi2dLYERERDnCzg4IDZW/T50K3L9v3niIiHIZowc9161bF6dOncqOWIiIiMyjWzegTh3g33+B8ePNHQ0RUa5i9BiGwYMHY9SoUbh79y5q1qwJR0dHjf1Vq1Y1WXBEREQ5wsICmDEDaNgQWLwY+OILoHJlc0dFRJQrKIQQwpgHWFhoN0ooFAoIIaBQKNQrP+dWiYmJcHV1xbNnz+Di4mLucIiI8pycrEdzvM7u1AlYvx4ICQF27sz+5yMiykamqkONbmG4ceNGpp+MiIgoV5s2DdiyBfjrL3kLCTF3REREZmd0wlCyZMnsiIOIiMj8ypSR3ZFmzJCLuTVtClgZ/VFJRJSvGD3oGQCuXbuGL774As2aNUPz5s3x5Zdf4tq1a6aOjYiIKOd99x1QuDBw/jywZIm5oyEiMjujE4a//voL/v7+OHbsGKpWrYrKlSvj6NGjqFSpEnbv3p0dMRIREeUcNzfghx/k799/DyQlmTceIiIzM3rQc40aNRASEoKpU6dqbP/666+xa9cunDx50qQBmhoHPRMRZU2+HvSs8vo1UKkScPUq8O23wKRJOffcREQmYqo61OgWhosXL+pc1XnAgAG4cOFCpgMhIiLKNWxsgJ9+kr//8gtw54554yEiMiOjE4ZixYohNjZWa3tsbCxXfSYiovyjfXugUSPg1Svgm2/MHQ0RkdkYPfXDp59+ioEDB+L69eto0KABFAoFDhw4gGnTpmHUqFHZESMREVHOUyhk60Lt2sCKFcCwYUCtWuaOiogoxxmdMHz//fdwdnbGL7/8gnHjxgEAvLy8MGHCBHz55ZcmD5CIiMhsatUCevWSCcOoUUBUlEwkiIgKEKMHPaeV9P8zRzg7O5ssoOzGQc9ERFlTIAY9p3XnDlCunOyatHEj0KGDeeIgIjKS2QY937hxA1euXAEgEwVVsnDlyhXcvHnT6ADmzZsHPz8/2NnZoWbNmoiOjk63fHJyMr799luULFkStra2KFOmDJZwnmwiIsouvr6ydQEAxoyRMygRERUgRicM/fr1w6FDh7S2Hz16FP369TPqWGvXrsXw4cPx7bff4tSpUwgMDETLli1x+/ZtvY/p0qUL9u7di99//x2XL1/G6tWrUaFCBWNfBhERkeHGjgXc3eU0q/PnmzsaIqIcZXSXJBcXF5w8eRLvvfeexvarV6+iVq1aePr0qcHHqlu3Lt5//33MT1P5VqxYER06dEBoaKhW+Z07d6Jbt264fv06ChcubEzYarmieZuIKA8rcF2SVBYtAj77TK4CffWqXOCNiCgXM1uXJIVCoR67kNazZ8+gVCoNPs7r169x4sQJtGjRQmN7ixYtdLZgAMCWLVtQq1YtTJ8+Hd7e3ihXrhxGjx6Nly9f6n2e5ORkJCYmatyIiCh3ytV19oABcjG3x4+5kBsRFShGJwyBgYEIDQ3VSA6USiVCQ0PRsGFDg4/z8OFDKJVKuLu7a2x3d3fH/fv3dT7m+vXrOHDgAM6dO4eNGzdi5syZWLduHYYMGaL3eUJDQ+Hq6qq++fr6GhwjERHlrFxdZ1tZyWlWAWD2bNnKQERUABjdJenChQto1KgRChUqhMDAQABAdHQ0EhMTERERgcqVKxt0nLi4OHh7e+PQoUOoX7++evvkyZMRFhaGS5cuaT2mRYsWiI6Oxv379+Hq6goA2LBhAzp16oTnz5/D3t5e6zHJyclITk5W309MTISvr2/uaN4mIsqDsrObUJ6osz/4APjrL+Cjj4B168wdDRGRXmbrkuTv748zZ86gS5cuSEhIQFJSEvr06YNLly4ZnCwAQNGiRWFpaanVmpCQkKDV6qDi6ekJb29vdbIAyDEPQgjcvXtX52NsbW3h4uKicSMiotwpT9TZP/8MWFgA69cDBw6YOxoiomxn9MJtgFyobcqUKVl6YhsbG9SsWRO7d+/Ghx9+qN6+e/dutG/fXudjAgICEB4ejn///RdOTk4AgL///hsWFhbw8fHJUjxEREQGqVwZ+Phj4Lff5HSrhw/LBIKIKJ8yuIZ7/Pix1lX88+fPo3///ujSpQtWrVpl9JOPHDkSixcvxpIlS3Dx4kWMGDECt2/fxueffw4AGDduHPr06aMu36NHDxQpUgT9+/fHhQsXsH//fowZMwYDBgzQ2R2JiIgoW0ycCDg5AceOAWvWmDsaIqJsZXDCMGTIEMyYMUN9PyEhAYGBgYiJiUFycjL69euHsLAwo568a9eumDlzJiZOnIjq1atj//792L59O0qWLAkAiI+P11iTwcnJCbt378bTp09Rq1Yt9OzZE23btsWvv/5q1PMSERFliYcH8PXX8vdx44B0ZusjIsrrDB707Ofnh6VLl6Jx48YAgJ9//hkLFizApUuXYGVlhZ9//hnr1q3DkSNHsjPeLMtVc3oTEeVBBXYdhne9eAGULw/cvQuEhr5NIIiIcokcH/R8//59+Pn5qe9HRETgww8/hJWVHAbRrl07XLlyJdOBEBER5SkODoBqPN/kycDGjcDq1UBUFGDEukRERLmdwQmDi4uLxirOx44dQ7169dT3FQqFxlR4RERE+V7PnkDp0sC//wIdOwI9egDBwUCpUsCGDeaOjojIJAxOGOrUqYNff/0VqampWLduHZKSktCkSRP1/r///jt3LbBDRESU3TZtAq5f195+7x7QqROTBiLKFwxOGP7zn/9g8+bNsLe3R9euXfHVV1/Bzc1NvX/NmjUICgrKliCJiIhyHaUSGDZM9z7V8MDhw9k9iYjyPIPXYahevTouXryIQ4cOwcPDA3Xr1tXY361bN/j7+5s8QCIiolwpOloOeNZHCODOHVnu/ycMISLKi4xauK1YsWJ6F1Vr3bq1SQIiIiLKE+LjDSs3bx5QpIhc8E2hyN6YiIiyAZemJCIiygxPT8PKhYcDVasCfn7AF18Au3YBnCSEiPIQJgxERESZERgI+PjobzVQKAA3N6BVK8DODrh1C5gzBwgJAYoVAzp3Bv74A3j4MGfjJiIyEhMGIiKizLC0BGbNkr+/mzSo7i9eDGzbJpOCzZuBTz4B3N2BpCRg3Tqgb195PzAQmD4duHjx7YDp9CiVcr0HrvtARDmACQMREVFmdewov/h7e2tu9/GR2zt2lPcdHYF27YDffgPi4oCjR4HvvgOqVQNSU4EDB4CxYwF/f6BcOWDkSCAyEnjzRvs5N2yQ6zwEB3PdByLKEQohDLmUoenatWtYunQprl27hlmzZqF48eLYuXMnfH19UalSpeyI02RMtUQ2EVFBlZP1aJ6ps5VKORtSfLwc2xAYKFsgDHHrFrB1K/DnnzJJeP367b5ChYCWLYG2bYEPPpD7O3XSboVQtWikTVJMKSuvj4jMxlR1qNEJw759+9CyZUsEBARg//79uHjxIkqXLo3p06fj2LFjWLduXaaDyQl55sOHiCiXYsKQjZKSgN27gS1b3nZlUrGwAKyt9Q+YVihky8aNG6b9Mr9hg1xvIu0Usj4+sjtWdiQnRGQyZksY6tevj86dO2PkyJFwdnbG6dOnUbp0acTExKBDhw64d+9epoPJCQXuw4eIyMSYMOQQpVJ2XfrzT3k7f96wx33zDVCzJuDgANjby5/v3uztASsDZlbfsME8LRpEZBJmSxicnJxw9uxZ+Pn5aSQMN2/eRIUKFfDq1atMB5MTCvSHDxGRCTBhMJOZM4ERI0x3PBsb7STi3fvbtwMvXuh+fHa1aBCRyZiqDjVq4TYAKFSoEOLj4+Hn56ex/dSpU/B+d9AXERERmUb16oaVq1lTTuP64gXw8qX8mfam8vq1vD19mrl4uJI1UYFhdMLQo0cPjB07FuHh4VAoFEhNTcXBgwcxevRo9OnTJztiJCIiItW6D/fu6Z56VXXF/+hR/Vf8hQBevdJMIHQlFS9eAPv3A8uXZxzXTz/JlayrVMna6yOiXMvohGHy5Mno168fvL29IYSAv78/lEolevToge+++y47YiQiIiLVug+dOsnkIG3SoBpTMHNm+t2DFArZ1cjeXn7JT0/p0oYlDNu3y9v778t1JXr0AIoWzfhxRJRnZGpaVQC4fv06Tp48idTUVNSoUQNly5Y1dWzZgv1hiYiyhmMYzEzXrEW+vjJZMOUAZKVSru+QXotGkSKy5WPr1rdrRlhbA61bA/36yVWura1NFxMRGcVsg57zOn74EBFlDROGXCCn1kVQzZIE6G7RUM2S9OiRXHV62TLgxIm35YoWBXr2lMmDoWMwiMhkTFWHGr3Sc6dOnTB16lSt7T/99BM6d+6c6UCIiIjIQJaWcqBx9+7yZ3bNUmToStZFigBDhwLHjwNnzwKjRwMeHnIdiVmzgBo15KrW//0vkJCQPbESUbYxuoWhWLFiiIiIQJV3BjedPXsWzZo1w4MHD0waoKnxahURUdawhaEAykyLRkoKsGuXbHXYvPntCtaWlrKrUr9+suuSra1pno+ItJhtWtV///0XNjY2Wtutra2RmJiY6UCIiIgol1K1aBjDykomBq1aAY8fA2vXykHUaRejK1xYDpLu21dOB6tQcGVpolzI6C5JlStXxtq1a7W2r1mzBv7+/iYJioiIiPKRwoWBQYOAI0eAixeBr78GvLxkIjFnDlC7tpyWtU8fOWYibbIAyIHXnTrJZIKIcpzRXZK2bNmCjz76CD169ECTJk0AAHv37sXq1asRHh6ODh06ZEecJsPmbSKirGGXJDIJpRLYu1d2Wdq4Ua4PkR6uLE1kNLMNem7Xrh02bdqEq1evYvDgwRg1ahTu3r2LPXv25PpkgYiIiHIJS0ugRQtg1Srg/n1g1Kj0y6ddWZqIcpTRYxgAoHXr1mjdurWpYyEiIqKCyNVVjmEwRHx89sZCRFoylTAAwOvXr5GQkIDU1FSN7SVKlMhyUERERFTAeHqathwRmYzRCcOVK1cwYMAAHDp0SGO7EAIKhQJKpdJkwREREVEBERgoxyjoW1kaANzcZDkiylFGJwz9+vWDlZUVtm7dCk9PTyhUqz0SERERZZalpZw6tVMnOcBZV9Lw5Anw44/yxu8fRDnG6IQhNjYWJ06cQIUKFbIjHiIiIiqoVCtLv7sOg68vUL8+8L//Af/5DxAXByxYINd6IKJsZ/R/mr+/Px4+fJgdsRAREVFB17Ej0L697pWemzaV6zn8/jvw4IFcDM7BwdwRE+V7Rk+rOm3aNHz11VeIiorCo0ePkJiYqHEjIiIiyhLVytLdu8ufqnUXBg6Ui7fZ2QFbt8oEghcxibKd0S0MzZo1AwA0bdpUYzsHPRMREVG2a99eLvjWtq1cObphQ2DnTqBUKXNHRpRvGZ0wREZGZkccRERERIZp0AA4cAD44APg8mU5vmHHDqB6dXNHRpQvGZ0wBAUFZUccRERERIarWBE4dAho2RI4exZo1AjYtAlo0sTckRHlO0aPYQCA6Oho9OrVCw0aNMC9e/cAAGFhYThw4IBJgyMiIiLSy9sb2L9fjnNISpItDmvWmDsqonzH6IRh/fr1CAkJgb29PU6ePInk5GQAQFJSEqZMmWLyAImIiIj0KlRIjmHo3Bl480YOlJ4509xREeUrRicMkyZNwoIFC/Dbb7/B2tpavb1BgwY4efKkSYMjIiIiypCtrWxZ+PJLeX/ECGDMGCA11bxxEeUTRicMly9fRqNGjbS2u7i44OnTp6aIiYiIiMg4FhayZWHaNHn/55+BPn2A16/NGhZRfmB0wuDp6YmrV69qbT9w4ABKly5tkqCIiIiIjKZQAF99Bfzxh1wFeuVKoE0bOb6BiDLN6IThs88+w7Bhw3D06FEoFArExcVh5cqVGD16NAYPHpwdMRIREREZrndvubCboyOwe7ccFH3/vrmjIsqzjJ5W9auvvsKzZ88QHByMV69eoVGjRrC1tcXo0aMxdOjQ7IiRiIiIyDghIUBUFNCqFXDypFy74a+/gLJlzR0ZoFQC0dFAfDzg6QkEBr5dzZooF1IIIYShhZVKJQ4cOIAqVarAzs4OFy5cQGpqKvz9/eHk5JSdcZpMYmIiXF1d8ezZM7i4uJg7HCKiPCcn61HW2ZRlV6/K6VavXQOKFgW2bQPq1DFfPBs2AMOGAXfvvt3m4wPMmgV07Gi+uPIqJl/pMlUdalSXJEtLS4SEhODZs2dwcHBArVq1UKdOnTyTLBAREVEB8957coG3mjWBhw+B4GC5KrQ5bNgAdOqkmSwAwL17cvuGDeaJK6/asAEoVUr+TXv0kD9LleJ5zAZGj2GoUqUKrl+/nh2xEBEREZle8eKye1JICPDiBdC2LbBsWc7GoFTKlgVdHTtU24YPl+UoY0y+cpTRCcPkyZMxevRobN26FfHx8UhMTNS4GWvevHnw8/ODnZ0datasiejoaL1lo6KioFAotG6XLl0y+nmJiIioAHFyAv78U061qlQC/fsDoaHyy7pSKROK1avlT1N8aX/1Cjh3Dti4UU712rat9pfbtIQA7twBdu3K+nPnd0y+cpzRg54/+OADAEC7du2gUCjU24UQUCgUUBrxx1m7di2GDx+OefPmISAgAAsXLkTLli1x4cIFlChRQu/jLl++rNEPq1ixYsa+DCIiIiporK1ly4KXFzB1KvDNN8D+/fKLfWbGFCiVwK1bwN9/a99u39b9hTYjrVsDVarIcRZ168qflSqxX35a0dGGJV/R0XKGLMoyoxOGyMhIkz35jBkz8PHHH+OTTz4BAMycORN//fUX5s+fj9DQUL2PK168OAoVKmSyOIiIiKiAUChky4KXl1wZeudO7TKqbi3r1gEffgg8eKA7Kbh2Lf2F4VxdgfLlgXLlZLKydGnG8QkBnDkjb4sXy22OjnIMhiqBqFtXJjVpLtymK78NDI6LM6zcoUNAo0ZyUT/KEqNmSTKl169fw8HBAeHh4fjwww/V24cNG4bY2Fjs27dP6zFRUVEIDg5GqVKl8OrVK/j7++O7775DcHCw3udJTk5GcnKy+n5iYiJ8fX054wYRUSZl58xFrLMpxyiVgLs78OiR/jLW1oCtLfDvv/rL2NrKqVrLldO+FS369ku9UikH5N67p7vlQaGQScDBg8CJE8DRo8CxY0BMjO6F5zw8NFshateWCcq78tOsTFevAitWAL/9ZnjSULQo0LQp0Lw50KwZULJk9saYy5iqvja6hQEAoqOjsXDhQly/fh3h4eHw9vZGWFgY/Pz80LBhQ4OO8fDhQyiVSri7u2tsd3d3x309i6t4enpi0aJFqFmzJpKTkxEWFoamTZsiKioKjRo10vmY0NBQ/Pjjj8a9QCIiMgvW2ZRjoqPTTxYA4M0bebOwkF/2dSUFPj6GXa23tJRf0jt1kslB2qRBlVTMnAn4+spbhw5ym1IJXL78NoE4elS2Pty/D2zZIm8qFSpotkJcuwZ066adoKRtQcntScOjR8DatUBYGHDkyNvt757Dd9nZyXP+8KF8/Nq1cnvZsjJ5aN5czqqkK8kiLUa3MKxfvx69e/dGz549ERYWhgsXLqB06dKYN28etm7diu3btxt0nLi4OHh7e+PQoUOoX7++evvkyZMRFhZm8EDmtm3bQqFQYEvaf5g0eLWKiMi02MJA+cLq1XIqzoz89BPwxReyJcEUdF3x9/WVyYKhX95fvpSL0akSiGPHgBs3jItD1aJx40bu65706pVcLyMsDNi+XSZtgEzcmjWTK3krFPInoDv5WrdODjQ/elSu9r17tzxPacfaWlrK5EqVQNStK1uV0pPHuneZrYVh0qRJWLBgAfr06YM1a9aotzdo0AATJ040+DhFixaFpaWlVmtCQkKCVqtDeurVq4cVK1bo3W9rawtbU/2TExFRtmKdTTnG09OwcrVqmS5ZAGRS0L591r502tsDAQHyppKQILsvqRKIgwfT70qlGhi8eXPuaGVITZUxh4UB4eHA06dv91WvLpOD7t01/2729rq7W6VNvho2lLcffwSePQMiI2XysGePHIdy+LC8TZwIODvLQdKq7ksVKmiOE8lP3buMZHQLg4ODAy5cuIBSpUrB2dkZp0+fRunSpXH9+nX4+/vj1atXBh+rbt26qFmzJubNm6fe5u/vj/bt26c76DmtTp064fHjx4iIiDCoPFcNJSLKGq70TPmCoWMKcuMVeEOsWgX07GlY2RIl3o6DqFNHDrB2ds7e+FQuX5ZJwsqVwM2bb7f7+Mj4e/UCKlfW//isXPG/dUsmDqoE4t0uaj4+MnFo3hxITgY+/lj7vZK2RSMXJg1ma2Hw9PTE1atXUapUKY3tBw4cQOnSpY061siRI9G7d2/UqlUL9evXx6JFi3D79m18/vnnAIBx48bh3r17+OOPPwDIWZRKlSqFSpUq4fXr11ixYgXWr1+P9evXG/syiIiIqCAzdExBXkwWADkLlKFu35a3devkfYUCqFhRJg+qW5UqgI1Nxscy5At8QgKwZo1MFI4ff7vd2Vn+PXr3BoKCDJvdyNIy81Onliwpk4CPP5YtHLGxb7svHTggWxKWLUt/kT8h5PkaPly2HOXV90sGjE4YPvvsMwwbNgxLliyBQqFAXFwcDh8+jNGjR+OHH34w6lhdu3bFo0ePMHHiRMTHx6Ny5crYvn07Sv7/CPb4+Hjcvn1bXf7169cYPXo07t27B3t7e1SqVAnbtm1Dq1atjH0ZREREVNB17Ci/JGfUrSUvCgyUryOjFpTTp+Xt2LG3szLdvg1cuCBvqi/Ltraya1DaloiyZTW/1KfXZadlS9n9acUKOZWtaiyBpSXwwQcySWjbFnBwyK4zkj4LC+D99+Vt7Fg5TiQ6WiYPGzfKAeT6FIB1HzI1req3336L//73v+ruR7a2thg9ejT+85//mDxAU2PzNhFR1rBLEuU7eWwgq8E2bJBX7AH9A4N1JUX378vEISbmbSLx5Il2OVdXmTzUri3P4U8/6Z+5yM5ODmZWqV1bdjfq1g0oXjxzry+nGDpAftUqOc4iFzFVHZrpdRhevHiBCxcuIDU1Ff7+/nBycsp0EDmJHz6UH+XXzzrKnZgwEOUhppiVSQjg+nXNVogTJzQTAEOULClbEnr2lAOK84qoKDkFa0YiI3NdC4PZE4a8ih8+lN8U4EkbyEyYMBDlMdlxVenNG+D8eZk8bNokpz/NyN69QJMmWXtec8hogDwgk7BcOEDebIOenz9/jqlTp2Lv3r1ISEhAamqqxv7r169nOhgiMo6qtTkvr8lDRETZLCsDg/WxtpZjGqpXB5ycDEsYHjwwbQw5Jb0B8iq9e+e6ZMGUjE4YPvnkE+zbtw+9e/eGp6cnFGnnpyWiHKNUypYFXfVWAZm0gYiIcgND17QwtFxupG+AvJOTXO9i3jxgwACgTBnzxZiNjE4YduzYgW3btiEg7WIhRJRjhHg701vaOktXuTt3gNmz5YWPIkVyLEQiIipIDJ2RKTAw52MzJV2L7tWtK9dqOHRI7j982HwzPWUjoxMGNzc3FC5cODtiIco3TNldNC5OTlN9/LgcY3b8uJzC2lAjRsibp6dc+6ZKlbc//f0zV69xkDUREanl9zUt0tLVvSs8XE7HeuYMMHCgXF8in/XAMXrQ84oVK7B582YsX74cDnkwg+IAOspuWRmE/OCBZmJw/Lj8Uv4uKys52UR600KruLvr7zaqUMjW07RJROXKcmptKz2XEzjImjjomYh0MsWMTHnV/v1yQLdSCfz6K/DFF+aOCEAOz5JUo0YNjbEKV69ehRACpUqVgrW1tUbZkydPZjqYnMAPH8pO+gYh65ry+p9/ZGKQNjnQ1cXIwgKoVAmoVevtrWpVOd4svUkbVC3AN24AL17IySzOnQPOnn37859/dL8OGxu5yOe7iURMDNC5s2Gvj/IvJgxEpFdBboKeOVM26VtZySlWGzY0d0Q5mzD8+OOPBh9w/PjxmQ4mJ/DDh7KLata19MYVuLoCTZvKJOHWLe39CoX8ol6rFlCzpvxZvbr+bkOZXZNHJSFBM4E4d07enj/XXV7f5BCqfaoEpaB8NhRUTBiIiHQQQi7ctnYt4OEBnDxp9oHeXIchk/jhQ9nF0HVd0ipf/m1iUKsWUKOGnHDBGKZuAU5NBW7e1G6NuHhR7svI/PlA//6Ara3xz015AxMGIiI9/v0XqFdPNus3bAhERMguAWZi9oThxIkTuHjxIhQKBfz9/VGjRo1MB5GT+OFDpvb0qZwUYeFCYPPmjMv36AF8+qlMDlxdTRNDTrQAh4UBffoYVtbaWnabql1bJkK1a8sB1vrGRWSkILdw50ZMGIiI0nHlivzwS0wEvvxSDvIzE7Mt3JaQkIBu3bohKioKhQoVghACz549Q3BwMNasWYNixYplOhii3E4I2ZXowAHg4EH58/x5/d10dPn0U9Ovn5Mda/K8y9fXsHLOzkBS0tvxGSr29jJJUiUQtWoB5crJMRrp4SBrIiLKU8qWlVfZ2reXA6Dr1pVXC/Mwo1sYunbtimvXriEsLAwVK1YEAFy4cAF9+/bFe++9h9WrV2dLoKbCq1UFU2avUKekyFnS0iYIcXHa5d57D2jQAPjzT9nikNEg5Lx4dVw1RiOjQdbXr8v1H44fl4OkVQO6k5K0H+PsLLtkpW2JKFXq7RgMYwaRmxJbNNLHFgYiIgN89x0webK8YnbkiGx6z2Fm65Lk6uqKPXv2oHbt2hrbjx07hhYtWuDp06eZDiYn8MOn4DHmCnVSEnD06NsE4cgR2R0xLSsr+SU3IODtzd397XNlZRBybpfZ15eaKltoVQlETAxw6hTw8qV22SJFZPLw/vvAokXAo0e6Y8muBIwtGhljwkBEZAClEmjdGvjrLzmH+fHjQKFCORqC2RIGZ2dnREdHo3r16hrbT506haCgICQmJmY6mJzAD5+CJaMr1AsXynEEqgQhNlZ7YK+rq2w9CAiQ45dq105/sbP8Pg21qV5fSgpw4YJmS8Tp08CbN8bFM3OmnPrazQ0oXFheyMnsejnmatHIa5gwEBEZ6NEjeRXs5k2gTRs52DGjvrgmZLaEoX379nj69ClWr14NLy8vAMC9e/fQs2dPuLm5YePGjZkOJifww6fgMGSaU11KlpSJgar1oFIl469g5/cuLdn1+pKT5YxMMTFyVrp9+4w/ho3N2+TBzc3w311c5KxV+t4veb1LmSkxYSAiMsLJk/ILxatXwI8/Aj/8kGNPbbaE4c6dO2jfvj3OnTsHX19fKBQK3L59G1WqVMHmzZvh4+OT6WByAj98Cg5Dpzl97z2gZcu3CUIufwsXGIb+/Xx8ZKLx5IlstchukZHZP8A8t2PCQERkpGXL5JzjCgWwbZv84pEDzDZLkq+vL06ePIndu3fj0qVLEELA398fzZo1y3QQRKaUmiqnOZ02zbDyEyfKdVYodwkMlMmAIStZW1rKMv/+KxOHJ0+Ax48N/13fQHVddC24R0RElK5+/eQgyQUL5IxJJ04ApUubOyqDceE2yhdSU+UA5fBwebt3z/DH8opx7pVTg8hTU+UFn3btMi5rZycvDLVrJ8eyFcSZpNnCQESUCcnJQFCQTByqVQMOHUp/UKQJmKoONXjURUREBPz9/XUOan727BkqVaqE6OjoTAdC5qVUyi4gq1fLn0qluSPKmKolYcQIOe4gIEAOgL13T07X2aOHnHFH3wBYhUIO1g0MzNGwyQgdO8qkwNtbc7uPj2kHIFtYAK1ayeOmN2Da0lJ2Qd24UbYsu7vL8S7TpwOXLhm3HgcRERUwtrbyw6t4cTnLx+ef55kPDoNbGNq1a4fg4GCMGDFC5/5ff/0VkZGRHPScB+WlaSRTU2VirmpJSBuzs7NcI6VzZ6BFC3klOL9Pc1pQ5NQg8ozeL+HhgJ8fsGWLvJ06pfn4smWBtm1l60NAQOZXts7t2MJARJQFUVFAs2byw23OHGDIkGx7qhwf9FyyZEns3LlTvVjbuy5duoQWLVrg9u3bmQ4mJ/DDR1NeWBhLCM0k4c6dt/ucnN4mCSEhMkl4V36f5pRMy5j3y507wNatMnmIiABev367z81Ndllq106+N/VVN3lxRi0mDEREWfTLL8Do0fLK0r59cv72bJDjCYOdnR3OnTuH9957T+f+q1evokqVKnipayWmXIQfPm9lNO2oORfGEgI4dgz43/9k0pI2D3Vykl/CunTRnyS8Ky9+KSPzycz7JSkJ2LVLJg9bt8oB1SrW1nLGp3btZAtEiRJye15q3UuLCQMRURYJAXTrJr/oeHnJQdAeHiZ/mhxPGMqUKYOff/4ZH374oc79GzZswOjRo3H9+vVMB5MT+OHzlqHTVlasCJQrJ7vcubvr/unmZtg6JBm1aEyZAjx8KFsSdCUJqpYEe3uDXyZRjktJkeNrtmyRa/RcuaK5v3p12X0pPFz7sXmhuxwTBiIiE/j3X6BuXbmKaaNGwJ498gqTCeV4wvDFF18gKioKMTExsHvnku7Lly9Rp04dBAcH49dff810MDmBHz5vLV4MfPqpaY5lZSVni9GXULi7ywHI7dvLq7aGcHR8myR88AGTBMq7Ll9+O+7h0CHt1cTfldsXiWPCQERkIpcvA7Vry2bqESOAGTNMevgcTxgePHiA999/H5aWlhg6dCjKly8PhUKBixcvYu7cuVAqlTh58iTc3d0zHUxO4IcPEBcn+2PPmQMY0oNs4kSZDCQkAA8eaP98+tS08QUHA0OHyqkrmSRQfvPPP7LrqiHrhHz3nVwjpFy53DWAmgkDEZEJbdz4tkl59WrZVclEzLLS861btzBo0CD89ddfUD1MoVAgJCQE8+bNQ6lSpTIdSE4pyB8+ly8DP/0EhIW9HZxpZaV/dVxDr3K+fi2TB1UCoSupSEgAbt4Enj3LOM5Vq7iQGuVvq1fLaX8NZWcHVKokuzJVqyZvVasChQoZ/9ymGM/DhIGIyMTGjQOmTpXrMhw9ClSubJLDmmWl55IlS2L79u148uQJrl69CiEEypYtCzc3t0wHQNnv2DF5NXPjxrdjBwIDgbFj5ZzynTvLbbqmkZw5M+MvEzY2MrHw8Um/nKFjJjw9My5DlJcZ+h7395crSz9/LsfDnTihub9kSZk8pE0k/Pz0jyfKq4OsiYjyvUmTgOPH5TiGjh2BmBjA1dXcUalxped8Sgjgr79kohAV9XZ7u3YyUUg7e1dOTTuqmpXp3j3d65Tk9n7bRKZizP+CQgFcvw7Exsp1flQ3fTNYOznJ1oe0iUTlyrI+MNUUymxhICLKBg8fAjVrygq+XTt5pdeQGWXSYZYuSflBfv/wSUmRM69Mmya/VACy21GvXsCYMfKKpS65ZWGs3DwzDJEpZfV/4ckT4MwZ+X+uSibOnweSk3WXN0X3QxUmDERE2eT4caBhQ1mZT5oEfPttlg7HhCGT8uuHz4sXwNKlcjDljRtym6MjMHCgHHTv62ve+NLiQmpEkqn/F1JS5FiltC0Rp08D9+8b9vjISKBx44zLMWEgIspGv/8OfPKJvJqzdasc15DJK7pMGDIpv334PH4MzJsH/PqrnH0FAIoWBb78Uq40XriweePThwupEUk58b+wYAEwaFDG5QydcIAJAxFRNhs4EPjtN5k0pP2qbuTAM7MMeqacld4Xibt35VS9ixbJAZGA7BM9ejTQv79MRnMzS0vDrmQS5Xc58b9QoYJh5TjhABFRLtG0qUwY3r2uf++e7M+aw324mTAYISeviuubzWTMGODUKWDlSuDNG7m9alU5kLlLl9w1VzsR5Q6BgbL+yGiQdWBgzsdGRETvUCrlFWBdhJCV9vDhcjXcHOqewa+XBsrJ6QhVgyHf/WC/e1fGoBIUJBOFDz54O1CSiOhdlpayrurUSbt125gplImIKAdER2t+4XyXEMCdO7JcDnXXyNpcTQWE6gv8u387VavQhg1Zf46UFLkqeHy8XOU4vZEl9vbAwYNyutSWLZksEFHGOnaULdje3prbfXw4OxkRUa4SH2/acibAFoYMKJXyqr6uL/CqbR9/DFy9KmfAevkyczd90x3q8vLl25WaiYgM1bGjbMHmhANERLmYoQPKcnDgGROGDGTUKgQAT5/KrkE5KQeTSiLKRzjhABFRLpcLB54xYciAoV/MAwLkomj29obdHBx0bz98WA6MzwhnMyEiIiLKh3LhwDMmDBkw9Iv5pEmmuWoXFJTrkkoiIiIiykmqgWe6Ztwxw0q3TBgykNOtQrkwqSQiIiKinJaLBp5xlqQMqL7AA9qzEWXXF3jOZkJERERE6oFn3bvLn2a6YsyEwQDm+ALfsSNw8yYQGQmsWiV/3rjBZIGIiIiIcha7JBnIHK1CnM2EiIiIiMytwCUM4v8HBSQmJmbq8e+///b3589NERERUd6iqj9FeitMmkhW62wiooLMVPV1gUsYkpKSAAC+vr5mjoSIKG9LSkqCq6trtj8HwDqbiCgrslpfK0ROXCLKRVJTUxEXFwdnZ2co3h3FnI8kJibC19cXd+7cgYuLi7nDyRV4TnTjedGN50Wb6pzcvn0bCoUCXl5esLDI3qFwBaHO5ntNN54X3XhetPGc6KY6LxcuXED58uWzVF8XuBYGCwsL+Pj4mDuMHOPi4sJ/nnfwnOjG86Ibz4s2V1fXHDsnBanO5ntNN54X3XhetPGc6Obt7Z3lizucJYmIiIiIiPRiwkBERERERHoxYcinbG1tMX78eNja2po7lFyD50Q3nhfdeF608ZxkD55X3XhedON50cZzopspz0uBG/RMRERERESGYwsDERERERHpxYSBiIiIiIj0YsJARERERER6MWHIw/bv34+2bdvCy8sLCoUCmzZt0tgvhMCECRPg5eUFe3t7NG7cGOfPnzdPsDkkNDQUtWvXhrOzM4oXL44OHTrg8uXLGmUK4nmZP38+qlatqp6jun79+tixY4d6f0E8J+8KDQ2FQqHA8OHD1dsK4nmZMGECFAqFxs3Dw0O9vyCeE1Nhna2NdbZurLMzxjpbyqk6mwlDHvb8+XNUq1YNc+bM0bl/+vTpmDFjBubMmYOYmBh4eHigefPmSEpKyuFIc86+ffswZMgQHDlyBLt370ZKSgpatGiB58+fq8sUxPPi4+ODqVOn4vjx4zh+/DiaNGmC9u3bqyuNgnhO0oqJicGiRYtQtWpVje0F9bxUqlQJ8fHx6tvZs2fV+wrqOTEF1tnaWGfrxjo7fayzNeVInS0oXwAgNm7cqL6fmpoqPDw8xNSpU9XbXr16JVxdXcWCBQvMEKF5JCQkCABi3759Qgiel7Tc3NzE4sWLC/w5SUpKEmXLlhW7d+8WQUFBYtiwYUKIgvteGT9+vKhWrZrOfQX1nGQH1tm6sc7Wj3W2xDpbU07V2WxhyKdu3LiB+/fvo0WLFupttra2CAoKwqFDh8wYWc569uwZAKBw4cIAeF4AQKlUYs2aNXj+/Dnq169f4M/JkCFD0Lp1azRr1kxje0E+L1euXIGXlxf8/PzQrVs3XL9+HUDBPifZjedWYp2tjXW2JtbZ2nKizrYyacSUa9y/fx8A4O7urrHd3d0dt27dMkdIOU4IgZEjR6Jhw4aoXLkygIJ9Xs6ePYv69evj1atXcHJywsaNG+Hv76+uNAriOVmzZg1OnjyJmJgYrX0F9b1St25d/PHHHyhXrhwePHiASZMmoUGDBjh//nyBPSc5geeWdfa7WGdrY52tLafqbCYM+ZxCodC4L4TQ2pZfDR06FGfOnMGBAwe09hXE81K+fHnExsbi6dOnWL9+Pfr27Yt9+/ap9xe0c3Lnzh0MGzYMu3btgp2dnd5yBe28tGzZUv17lSpVUL9+fZQpUwbLly9HvXr1ABS8c5KTCvK5ZZ2tiXW2JtbZuuVUnc0uSfmUaoS8KrtUSUhI0Mo086MvvvgCW7ZsQWRkJHx8fNTbC/J5sbGxwXvvvYdatWohNDQU1apVw6xZswrsOTlx4gQSEhJQs2ZNWFlZwcrKCvv27cOvv/4KKysr9WsvaOflXY6OjqhSpQquXLlSYN8rOaGgn1vW2dpYZ2tinW2Y7KqzmTDkU35+fvDw8MDu3bvV216/fo19+/ahQYMGZowsewkhMHToUGzYsAERERHw8/PT2F9Qz4suQggkJycX2HPStGlTnD17FrGxsepbrVq10LNnT8TGxqJ06dIF8ry8Kzk5GRcvXoSnp2eBfa/khIJ6bllnG451NutsQ2RbnW3UEGnKVZKSksSpU6fEqVOnBAAxY8YMcerUKXHr1i0hhBBTp04Vrq6uYsOGDeLs2bOie/fuwtPTUyQmJpo58uwzaNAg4erqKqKiokR8fLz69uLFC3WZgnhexo0bJ/bv3y9u3Lghzpw5I7755hthYWEhdu3aJYQomOdEl7QzbghRMM/LqFGjRFRUlLh+/bo4cuSIaNOmjXB2dhY3b94UQhTMc2IqrLO1sc7WjXW2YVhn51ydzYQhD4uMjBQAtG59+/YVQsjptMaPHy88PDyEra2taNSokTh79qx5g85mus4HALF06VJ1mYJ4XgYMGCBKliwpbGxsRLFixUTTpk3VHzxCFMxzosu7Hz4F8bx07dpVeHp6Cmtra+Hl5SU6duwozp8/r95fEM+JqbDO1sY6WzfW2YZhnZ1zdbZCCCEy2epBRERERET5HMcwEBERERGRXkwYiIiIiIhILyYMRERERESkFxMGIiIiIiLSiwkDERERERHpxYSBiIiIiIj0YsJARERERER6MWEgIiIiIiK9mDAQGejmzZtQKBSIjY01dyhqly5dQr169WBnZ4fq1avn+PPnxnNCRATkzvqJdTblVUwYKM/o168fFAoFpk6dqrF906ZNUCgUZorKvMaPHw9HR0dcvnwZe/fu1dqvUCjSvfXr1y9Lz+/r64v4+HhUrlw5S8chovyHdbY21tmUVzFhoDzFzs4O06ZNw5MnT8wdism8fv0604+9du0aGjZsiJIlS6JIkSJa++Pj49W3mTNnwsXFRWPbrFmzshI6LC0t4eHhASsrqywdh4jyJ9bZmlhnU17FhIHylGbNmsHDwwOhoaF6y0yYMEGrqXfmzJkoVaqU+n6/fv3QoUMHTJkyBe7u7ihUqBB+/PFHpKSkYMyYMShcuDB8fHywZMkSreNfunQJDRo0gJ2dHSpVqoSoqCiN/RcuXECrVq3g5OQEd3d39O7dGw8fPlTvb9y4MYYOHYqRI0eiaNGiaN68uc7XkZqaiokTJ8LHxwe2traoXr06du7cqd6vUChw4sQJTJw4EQqFAhMmTNA6hoeHh/rm6uoKhUKhsW3VqlUoU6YMbGxsUL58eYSFhWk8XqFQYP78+WjZsiXs7e3h5+eH8PBw9X5dzdvnz59H69at4eLiAmdnZwQGBuLatWsAgKioKNSpUweOjo4oVKgQAgICcOvWLZ2vn4jyPtbZrLMpf2DCQHmKpaUlpkyZgtmzZ+Pu3btZOlZERATi4uKwf/9+zJgxAxMmTECbNm3g5uaGo0eP4vPPP8fnn3+OO3fuaDxuzJgxGDVqFE6dOoUGDRqgXbt2ePToEQB5dSgoKAjVq1fH8ePHsXPnTjx48ABdunTROMby5cthZWWFgwcPYuHChTrjmzVrFn755Rf8/PPPOHPmDEJCQtCuXTtcuXJF/VyVKlXCqFGjEB8fj9GjRxv1+jdu3Ihhw4Zh1KhROHfuHD777DP0798fkZGRGuW+//57fPTRRzh9+jR69eqF7t274+LFizqPee/ePTRq1Ah2dnaIiIjAiRMnMGDAAKSkpCAlJQUdOnRAUFAQzpw5g8OHD2PgwIEFtmsCUUHAOpt1NuUTgiiP6Nu3r2jfvr0QQoh69eqJAQMGCCGE2Lhxo0j7Vh4/fryoVq2axmP/+9//ipIlS2ocq2TJkkKpVKq3lS9fXgQGBqrvp6SkCEdHR7F69WohhBA3btwQAMTUqVPVZd68eSN8fHzEtGnThBBCfP/996JFixYaz33nzh0BQFy+fFkIIURQUJCoXr16hq/Xy8tLTJ48WWNb7dq1xeDBg9X3q1WrJsaPH5/hsYQQYunSpcLV1VV9v0GDBuLTTz/VKNO5c2fRqlUr9X0A4vPPP9coU7duXTFo0CAhxNtzcurUKSGEEOPGjRN+fn7i9evXWs//6NEjAUBERUUZFC8R5W2ss1lnU/7BFgbKk6ZNm4bly5fjwoULmT5GpUqVYGHx9l/A3d0dVapUUd+3tLREkSJFkJCQoPG4+vXrq3+3srJCrVq11FdvTpw4gcjISDg5OalvFSpUAAB1Ey8A1KpVK93YEhMTERcXh4CAAI3tAQEBeq8UGevixYsGHT/t61Xd1xdDbGwsAgMDYW1trbWvcOHC6NevH0JCQtC2bVvMmjUL8fHxWXwVRJQXsM7OOtbZZE5MGChPatSoEUJCQvDNN99o7bOwsIAQQmPbmzdvtMq9W0EqFAqd21JTUzOMR9VEm5qairZt2yI2NlbjduXKFTRq1Ehd3tHRMcNjpj2uihDCpM3BmT2+vjL29vbpPm7p0qU4fPgwGjRogLVr16JcuXI4cuSI4QETUZ7EOts0WGeTuTBhoDxr6tSp+PPPP3Ho0CGN7cWKFcP9+/c1PoBMOed02soyJSUFJ06cUF+Rev/993H+/HmUKlUK7733nsbN0A8cAHBxcYGXlxcOHDigsf3QoUOoWLGiSV5HxYoVDTr+ux8OR44cUb/ed1WtWhXR0dE6P+xVatSogXHjxuHQoUOoXLkyVq1alclXQER5CevsrGGdTebEhIHyrCpVqqBnz56YPXu2xvbGjRvjn3/+wfTp03Ht2jXMnTsXO3bsMNnzzp07Fxs3bsSlS5cwZMgQPHnyBAMGDAAADBkyBI8fP0b37t1x7NgxXL9+Hbt27cKAAQOgVCqNep4xY8Zg2rRpWLt2LS5fvoyvv/4asbGxGDZsmElex5gxY7Bs2TIsWLAAV65cwYwZM7BhwwatgXjh4eFYsmQJ/v77b4wfPx7Hjh3D0KFDdR5z6NChSExMRLdu3XD8+HFcuXIFYWFhuHz5Mm7cuIFx48bh8OHDuHXrFnbt2oW///7bZB+mRJS7sc7OGtbZZE5MGChP+89//qPVlF2xYkXMmzcPc+fORbVq1XDs2DGjZ6NIz9SpUzFt2jRUq1YN0dHR2Lx5M4oWLQoA8PLywsGDB6FUKhESEoLKlStj2LBhcHV11eh7a4gvv/wSo0aNwqhRo1ClShXs3LkTW7ZsQdmyZU3yOjp06IBZs2bhp59+QqVKlbBw4UIsXboUjRs31ij3448/Ys2aNahatSqWL1+OlStXwt/fX+cxixQpgoiICPz7778ICgpCzZo18dtvv8Ha2hoODg64dOkSPvroI5QrVw4DBw7E0KFD8dlnn5nk9RBR7sc6O/NYZ5M5KcS7/7lERP9PoVBg48aN6NChg7lDISKiDLDOpuzCFgYiIiIiItKLCQMREREREenFLklERERERKQXWxiIiIiIiEgvJgxERERERKQXEwYiIiIiItKLCQMREREREenFhIGIiIiIiPRiwkBERERERHoxYSAiIiIiIr2YMBARERERkV5MGIiIiIiISK//A/nYpwv0WJCRAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 800x300 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot coherence scores\n",
    "topic_numbers = list(range(4, 49, 4))\n",
    "\n",
    "# Create a figure with two subplots\n",
    "fig, axs = plt.subplots(1, 2, figsize=(8, 3), sharey=True)\n",
    "\n",
    "# Plot coherence values for each category\n",
    "axs[0].plot(topic_numbers, coherence_values_0, marker='o', linestyle='-', color='b')\n",
    "axs[0].set_title('June')\n",
    "axs[0].set_xlabel('Number of Topics')\n",
    "axs[0].set_ylabel('Coherence Score')\n",
    "\n",
    "axs[1].plot(topic_numbers, coherence_values_1, marker='o', linestyle='-', color='r')\n",
    "axs[1].set_title('Oct/Nov')\n",
    "axs[1].set_xlabel('Number of Topics')\n",
    "# axs[1].set_ylabel('Coherence Score')  # Shared y-axis\n",
    "\n",
    "plt.suptitle('Coherence Scores for LDA Models')\n",
    "plt.tight_layout(rect=[0, 0.03, 1, 0.95]) # Adjust layout to make room for the main title\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "75a2c5df-77c0-406b-bb11-345d606a80c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_entropy(model, corpus):\n",
    "    \"\"\"\n",
    "    Calculate the entropy of the topic distribution for a given LDA model and corpus.\n",
    "    \n",
    "    Entropy is a measure of the unpredictability or disorder within the topic distribution\n",
    "    for each document in the corpus, and this function calculates the average entropy across\n",
    "    all documents to provide a measure of the overall diversity of topics within the corpus.\n",
    "    \n",
    "    Parameters:\n",
    "    - model: An instance of LdaModel trained on the corpus.\n",
    "    - corpus: The corpus on which the model is trained, in the gensim bag-of-words format.\n",
    "    \n",
    "    Returns:\n",
    "    - avg_entropy: The average entropy of the topic distribution across all documents in the corpus.\n",
    "    \"\"\"\n",
    "    entropy_list = []\n",
    "    for document in corpus:\n",
    "        topic_distribution = model.get_document_topics(document, minimum_probability=0)\n",
    "        probabilities = np.array([probability for _, probability in topic_distribution])\n",
    "        entropy = -np.sum(probabilities * np.log(probabilities))\n",
    "        entropy_list.append(entropy)\n",
    "    \n",
    "    avg_entropy = np.mean(entropy_list)\n",
    "    return avg_entropy\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "001b725f-ca53-4307-af14-f20a6f92244a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.29249054, 0.45124906)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate the entropy of the topic distribution\n",
    "entropy_category_0 = calculate_entropy(model_list_0[6], corpus_0)\n",
    "entropy_category_1 = calculate_entropy(model_list_1[6], corpus_1)\n",
    "\n",
    "entropy_category_0, entropy_category_1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "0b0184c2-27a9-4436-a91f-3537929afe4b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0.19852467, 0.25641456),\n",
       " (0.30142248, 0.366717),\n",
       " (0.32976708, 0.439908),\n",
       " (0.35311833, 0.4630826),\n",
       " (0.33727166, 0.43244916),\n",
       " (0.2893856, 0.45624092),\n",
       " (0.29249132, 0.45122057),\n",
       " (0.29130095, 0.47043666),\n",
       " (0.30828246, 0.45271713),\n",
       " (0.3047557, 0.49212164),\n",
       " (0.25659448, 0.51710653),\n",
       " (0.32417554, 0.50778174)]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate entropies for all numbers of topics\n",
    "entropies = [(calculate_entropy(model_list_0[i], corpus_0),calculate_entropy(model_list_1[i], corpus_1)) for i in range(len(model_list_0))]\n",
    "entropies"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "catching_trolls",
   "language": "python",
   "name": "catching_trolls"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
